options(warn = -1)
suppressPackageStartupMessages({
library(plotly)
library(dplyr)
library(glmnet)
library(corrplot)
library(caret)
})
Nota: Nos aseguramos de tener el paquete
ElemStatLearn instalado. Si no:
install.packages(lib_URL)
lib_URL = 'https://cran.r-project.org/src/contrib/Archive/ElemStatLearn/ElemStatLearn_2015.6.26.2.tar.gz'
# install.packages(lib_URL)
# install.packages("plotly")
library(plotly)
library(ElemStatLearn)
library(dplyr)
Utilizamos head() para ver las primeras filas y hacernos
una idea de las variables disponibles
data("prostate", package = "ElemStatLearn")
head(prostate)
## lcavol lweight age lbph svi lcp gleason pgg45 lpsa
## 1 -0.5798185 2.769459 50 -1.386294 0 -1.386294 6 0 -0.4307829
## 2 -0.9942523 3.319626 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 3 -0.5108256 2.691243 74 -1.386294 0 -1.386294 7 20 -0.1625189
## 4 -1.2039728 3.282789 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 5 0.7514161 3.432373 62 -1.386294 0 -1.386294 6 0 0.3715636
## 6 -1.0498221 3.228826 50 -1.386294 0 -1.386294 6 0 0.7654678
## train
## 1 TRUE
## 2 TRUE
## 3 TRUE
## 4 TRUE
## 5 TRUE
## 6 TRUE
attach(prostate)
TRAIN.La columna train se utiliza en estudios previos para
indicar las filas de entrenamiento.
Dado que esta práctica no la requiere, eliminamos esta columna para simplificar el conjunto de datos.
prostate = prostate %>% select(-train)
head(prostate)
## lcavol lweight age lbph svi lcp gleason pgg45 lpsa
## 1 -0.5798185 2.769459 50 -1.386294 0 -1.386294 6 0 -0.4307829
## 2 -0.9942523 3.319626 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 3 -0.5108256 2.691243 74 -1.386294 0 -1.386294 7 20 -0.1625189
## 4 -1.2039728 3.282789 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 5 0.7514161 3.432373 62 -1.386294 0 -1.386294 6 0 0.3715636
## 6 -1.0498221 3.228826 50 -1.386294 0 -1.386294 6 0 0.7654678
Tras eliminar la columna train, contamos con un total de
9 variables.
Utilizamos la función sapply para aplicar una función,
en este caso class, a todas las columnas.
sapply(prostate, class)
## lcavol lweight age lbph svi lcp gleason pgg45
## "numeric" "numeric" "integer" "numeric" "integer" "numeric" "integer" "integer"
## lpsa
## "numeric"
Los tipos de las variables son del tipo numeric e
integer, por lo que todos datos son números. Dado que son
datos numéricos, podemos echar un vistazo a sus distribuciones.
par(mfrow=c(3,3))
for (column in names(prostate)) {
hist(prostate[[column]], main=paste("Histograma de", column), xlab = column)
}
Podemos observar que los valores de la variable svi se
concentran en el 0 y en el 1, indicador de que esta variable podría ser
un booleano. Veamos si estos son sus únicos valores.
unique(prostate$svi)
## [1] 0 1
0 y 1 son los únicos valores presentes, por lo que confirmaríamos que se trata de una variable booleana donde 1 indicaría que el cáncer ha invadido el vesículo seminal, y 0 indicaría lo contrario.
No existe ninguna variable con el identificador del paciente.
Averiguamos si existe algún valor nulo en todo el conjunto de datos. En caso de que lo haya, identificamos las columnas que los contienen.
any(is.na(prostate))
## [1] FALSE
Dado que no hay ningún valor nulo en nuestro set de datos, continuamos.
Tras revisar los histogramas realizados en el subapartado 2 podemos concluir que, dados los rangos de las variables y su disposición, estos no se encuentran normalizados.
El objetivo de poner algunas de las variables en escala logarítmica es el de facilitar su posterior análisis por varios motivos:
lcavol, lweight y
lpsa).Además, cabe añadir que en los campos de la biología y la química es común usar escalas logarítmicas cuando estas hacen más fácil la interpretación de los datos.
Para convertir las variables svi y gleason
a categóricas usamos la función factor.
prostate$svi <- factor(prostate$svi)
prostate$gleason <- factor(prostate$gleason)
A continuación, visualizamos las distribuciones usando gráficos de barras:
par(mfrow = c(1, 2)) # Ajustar disposición gráfica
barplot(table(prostate$svi), main = "Distribución de SVI", col = "lightblue")
barplot(table(prostate$gleason), main = "Distribución de Gleason", col = "lightgreen")
Se puede observar como para la mayoría de pacientes se indica que el cáncer no ha invadido la vesícula seminal, a la par que la mayoría tiene un índice de Gleason de 6 o 7. Esto puede indicar cierta relación, ya que el índice de Gleason mide la agresividad del cáncer, por lo que índices más bajos indican cánceres menos agresivos. Podríamos aventurarnos a pensar que la mayoría de cánceres agresivos no llegarían a infectar la vesícula seminal.
En el caso de la variable Edad queremos hacer una
agrupación cada 5 años, por lo que usaremos la función cut
para generar esas categorías:
max.age = max(prostate$age)
min.age = min(prostate$age)
# Añadimos +5 al maximo para que el rango tome tambien el ultimo grupo
prostate$age <- cut(prostate$age, breaks = seq(from=min.age, to=max.age+5, by=5), right=FALSE)
barplot(table(prostate$age), main = "Distribución de Edad", col = "salmon", width = 3)
En el caso de la edad nos encontramos una distribución que podría recordar a la normal, centrada en los 65 años, encontrándose el grueso de la distribución entre los 61 y 70 años. Cabe añadir que encontramos casos desde los 41 hasta los 79 años.
# Filtramos por gleason = 7
gleason.7 <- subset(prostate, gleason == 7)
gleason.7.props <- prop.table(table(gleason.7$svi))*100
Por lo tanto, un 66% de los pacientes con índice de Gleason igual a 7
presenta un svi de 0. Podemos representar esta distribución
mediante un gráfico de sectores:
fig <- plot_ly(
labels = names(gleason.7.props),
values = as.numeric(gleason.7.props),
type="pie",
textinfo = "label+percent",
hoverinfo = "none",
insidetextorientation = "radial",
height = 400
) %>%
layout(
title = "Distribución del índice de Gleason para pacientes con svi = 0",
showlegend = TRUE
)
fig
# Filtramos por gleason = 7
svi.0 <- subset(prostate, svi == 0)
svi.0.props <- prop.table(table(svi.0$gleason))*100
Por lo tanto, un 49% de los pacientes con un svi de 0
presenta un índice de Gleason igual a 7. Podemos visualizar la
distribución del índice de Gleason para pacientes con un
svi\(=0\) mediante un
gráfico de sectores:
fig <- plot_ly(
labels = names(svi.0.props),
values = as.numeric(svi.0.props),
type="pie",
textinfo = "label+percent",
hoverinfo = "none",
insidetextorientation = "radial" ,
height = 400
) %>%
layout(
title = "Distribución del índice de Gleason para pacientes con svi = 0",
showlegend = TRUE
)
fig
Vistos los resultados de estos dos últimos subapartados, podemos pensar que estas variables están relacionadas. Además, ya en el apartado 2 se vio como también podría existir una posible relación entre las variables y se le trató de dar una explicación lógica.
Para probar si realmente existe esta relación entre variables, exponemos los datos a un test de independencia chi cuadrado, cuya implementación en R es simple y rápida. En este test queremos probar que las variables están relacionadas, lo cual es nuestra hipótesis \(H_1\), por lo que la hipótesis nula \(H_0\) será que estas son independientes:
contingencia.svi_gleason <- table(prostate$gleason, prostate$svi)
chi_test <- chisq.test(contingencia.svi_gleason)
chi_test
##
## Pearson's Chi-squared test
##
## data: contingencia.svi_gleason
## X-squared = 15.918, df = 3, p-value = 0.001179
Al realizar el test, obtenemos un \(p\)-valor de 0.0012. Dado que los niveles de significación más comunes son \(0.05\) y \(0.01\), y nuestro \(p\)-valor es notablemente menor, podemos afirmar con seguridad que existe una relación entre estas variables.
Análisis de la dependencia lineal de la variable lpsa
con la variable lcavol.
library(dplyr)
library(ggplot2)
Ajustamos el modelo y vemos su summary():
modelo_lineal <- lm(lpsa ~ lcavol, data = prostate)
summary(modelo_lineal)
##
## Call:
## lm(formula = lpsa ~ lcavol, data = prostate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.67624 -0.41648 0.09859 0.50709 1.89672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50730 0.12194 12.36 <2e-16 ***
## lcavol 0.71932 0.06819 10.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7875 on 95 degrees of freedom
## Multiple R-squared: 0.5394, Adjusted R-squared: 0.5346
## F-statistic: 111.3 on 1 and 95 DF, p-value: < 2.2e-16
Ecuación del Modelo:
La ecuación del modelo de regresión lineal es:
\[ \text{lpsa} = 1.50730 + 0.71932 \cdot \text{lcavol} \]
Intercepto (\(\beta_0\)):
lpsa cuando lcavol es igual a 0.Pendiente (\(\beta_1\)):
lcavol es 0.71932.
Esto sugiere que, por cada unidad de incremento en lcavol,
se espera que lpsa aumente en aproximadamente
0.71932 unidades.Crearemos una gráfico de dispersión de lcavol vs
lpsa.
ggplot(prostate,
aes(x = lcavol, y = lpsa)) +
# Puntos de los datos
geom_point(color = "blue",
alpha = 0.4) +
# Recta de regresión
geom_smooth(method = "lm",
color = "red") +
labs(title = "Relación entre lcavol y lpsa",
x = "Log de Volumen de Cáncer (lcavol)",
y = "Log de Antígeno Prostático Específico (lpsa)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
Los intervalos de confianza al 95% para los coeficientes del modelo son importantes para entender la precisión de las estimaciones. A continuación se presentan los intervalos calculados:
confint(modelo_lineal, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 1.2652222 1.7493727
## lcavol 0.5839404 0.8547004
Con esto podemos estimar que con ausencia de un tumor medido
(lcavol), el psa se encontrará entre valores
de 1.26 y 1.74 ng/m, en su medida logarítmica.
Ademas, podemos estimar que por cada cm (10mm) de tumor aumentado, el PSA aumenta en unidades comprendidos entre 5.8 y 8.5 ng/m.
En el análisis de regresión lineal, se ha observado que tanto el
coeficiente del intercepto como el coeficiente de lcavol
presentan valores \(p\) muy pequeños
\((< 2e-16)\). Esto indica que
podemos rechazar la hipótesis nula \((H_0)\) de que los “verdaderos”
coeficientes son nulos. Específicamente, esto significa que existe una
relación estadísticamente significativa entre el volumen del cáncer
(lcavol) y el nivel de antígeno específico de la próstata
(lpsa).
El Error Estándar de los Residuos (RSE) es una medida de la cantidad de variabilidad de los residuos (la diferencia entre los valores observados y los valores predichos por el modelo) en un modelo de regresión lineal. Se utiliza para evaluar la precisión de las predicciones del modelo. Un RSE más bajo indica que el modelo tiene un mejor ajuste a los datos.
Matemáticamente, el RSE se define como:
\[ RSE = \sqrt{\frac{1}{n-2} \sum (y_i - \hat{y}_i)^2} \]
donde:
\(y_i\) son los valores observados,
\(\hat{y}_i\) son los valores predichos por el modelo,
\(n\) es el número de observaciones.
Si echamos un vistazo de nuevo a nuestro resultado:
summary(modelo_lineal)
##
## Call:
## lm(formula = lpsa ~ lcavol, data = prostate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.67624 -0.41648 0.09859 0.50709 1.89672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50730 0.12194 12.36 <2e-16 ***
## lcavol 0.71932 0.06819 10.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7875 on 95 degrees of freedom
## Multiple R-squared: 0.5394, Adjusted R-squared: 0.5346
## F-statistic: 111.3 on 1 and 95 DF, p-value: < 2.2e-16
Un RSE de 0.7875 indica que, en promedio, las predicciones del modelo
se desvían de los valores observados de la variable de respuesta
(lpsa) en aproximadamente 0.7875 unidades en la escala
logarítmica.
Indicando que por cada cm(10mm) de aumento en el lcavol
tenemos un desvía de 7.8 ng/m.
El coeficiente de determinación \(R^2\) es una métrica que indica la
proporción de la varianza en la variable dependiente lpsa
que es predecible a partir de la variable independiente
lcavol. En este modelo, el \(R^2\) es:
\(R^2 = 0.5394\)
El \(R^2\) ajustado, que en este caso es:
\(R^2\) Ajustado = 0.5346
psa y
cavol)El modelo de regresión lineal simple que hemos calculado es:
\[ \text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol} \]
\[ \text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol} \]
\[ ln(\text{psa}) = 1.50730 + 0.71932 \times ln(\text{cavol}) \]
\[ \text{PSA} = e^{ln(\text{psa})} = e^{(1.50730 + 0.71932 \times ln(\text{cavol}))} \]
\[ \text{PSA} = e^{1.50730 } + e^{(0.71932 \times ln(\text{cavol}))} \]
\[ \text{PSA} = e^{1.50730 } + (e^{ln(\text{cavol})})^{0.71932} \]
cavol: Para expresar
esto en términos de cavol (el volumen de cáncer en su forma
original), tienes que revertir el logaritmo:\[ \text{cavol} = e^{ln(\text{cavol})} \]
Ahora sustituimos lcavol por su expresión en términos de
cavol:
\[ \text{PSA} = e^{1.50730 } +\text{cavol}^{0.71932} \]
PSA
y cavol queda como:\[
\text{PSA} = C \cdot \text{cavol}^{b}
\] Donde \(C = e^{1.50730}\) y
\(b = 0.71932\). Esto significa que
PSA es proporcional a cavol elevado a un
exponente.
PSA cuando cavol es igual a 1 (o
la unidad de referencia).cavol), el PSA aumenta en
un porcentaje específico que depende del valor de
cavol.library(corrplot)
prostate$svi = as.numeric(prostate$svi)
prostate2 = prostate
prostate2$age_num <- as.numeric(prostate2$age)
prostate2$pgg45_num <- as.numeric(prostate2$pgg45)
prostate2$gleason_num <- as.numeric(prostate2$gleason)
cor_data = prostate2 %>% select(-age, -pgg45, -gleason)
cor_matrix = cor(cor_data)
corrplot(cor_matrix,
diag = FALSE,
method = "circle",
type = "upper",
tl.col = "black",
tl.srt = 25,
addCoef.col = "black")
La matriz de correlación muestra las correlaciones de Pearson entre las variables del conjunto de datos, con valores que van de -1 a 1.
Se ha realizado un análisis de correlación entre las variables del conjunto de datos. A continuación se presentan los resultados y las interpretaciones:
lcavol, lpsa)
lcavol) y el
nivel de antígeno específico de la próstata (PSA) sugiere que a medida
que aumenta el volumen del cáncer, también tienden a aumentar los
niveles de PSA.lbph, lweight):
lweight) y la hipertrofia benigna de
próstata (lbph) indica que hay una relación positiva entre
ambos. Esto sugiere que los hombres con mayor peso pueden tener mayor
riesgo de desarrollar condiciones benignas en la próstata.lcp:
lcavol
(0.675).lcp) y el
volumen del cáncer (lcavol) tienen trivial relación.svi:
lcp (0.68).svi) y el nivel de extensión del cáncer
dentro de la cápsula (lcp) es positiva y moderada (0.68).
Esto indica que a medida que aumenta la invasión a las vesículas
seminales, también podría aumentar la extensión del cáncer dentro de la
cápsula.age_num:
lbph (0.364).age_num) y la hipertrofia
prostática benigna (lbph) sugiere que a mayor edad, la
probabilidad de presentar hipertrofia benigna aumenta.pgg45_num, gleason_num):
pgg45_num) y la
puntuación total de Gleason (gleason_num) refuerza que un
mayor porcentaje de patrones más agresivos se asocia con una puntuación
de Gleason más alta.lpsa.modelo_lpsa = lm(lpsa ~ lcavol + lweight + lbph + lcp + svi,
data = prostate)
summary(modelo_lpsa)
##
## Call:
## lm(formula = lpsa ~ lcavol + lweight + lbph + lcp + svi, data = prostate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.84878 -0.38372 -0.00413 0.45189 1.55468
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.12471 0.73259 -1.535 0.12819
## lcavol 0.54790 0.08637 6.343 8.56e-09 ***
## lweight 0.53036 0.19769 2.683 0.00867 **
## lbph 0.07999 0.05643 1.418 0.15971
## lcp -0.03638 0.08088 -0.450 0.65391
## svi 0.75975 0.24122 3.150 0.00221 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7071 on 91 degrees of freedom
## Multiple R-squared: 0.6443, Adjusted R-squared: 0.6248
## F-statistic: 32.97 on 5 and 91 DF, p-value: < 2.2e-16
Siendo la siguiente su función:
\[ \text{lpsa} = -1.12471 + 0.54790 \cdot \text{lcavol} + 0.53036 \cdot \text{lweight} + 0.07999 \cdot \text{lbph} - 0.03638 \cdot \text{lcp} + 0.75975 \cdot \text{svi} + \epsilon \]
lcavol: 0.54790 (muy significativo,
p = 8.56e-09). Esto indica que, manteniendo las otras
variables constantes, un aumento de una unidad en el volumen del cáncer
(lcavol) se asocia con un aumento de 0.54790 unidades en
los niveles de PSA (lpsa).
lweight: 0.53036 (significativo,
p = 0.00867). Similar a lcavol, un aumento de
una unidad en el peso del paciente (lweight) se asocia con
un aumento en los niveles de PSA, sugiriendo que el peso del paciente
podría influir en el diagnóstico.
lbph: 0.07999 (no significativo,
p = 0.15971). La hipertrofia benigna de próstata
(lbph) no tiene un efecto significativo en los niveles de
PSA, ya que su valor de p es mayor que 0.05.
lcp: -0.03638 (no significativo,
p = 0.65391). El nivel de extensión del cáncer dentro de la
cápsula (lcp) no tiene un efecto significativo sobre los
niveles de PSA.
svi: 0.75975 (muy significativo,
p = 0.00221). La invasión seminal (svi) tiene
una asociación significativa con los niveles de PSA, indicando que a
medida que la invasión seminal aumenta, también aumentan los niveles de
PSA.
\(R^{2}\) ajustado: 0.6248. El modelo explica aproximadamente el 62.48% de la variación en los niveles de PSA. Esto sugiere una mejora en la capacidad predictiva con la adición de la variable svi, aunque aún hay un 37.52% de la variación que no se explica por las variables incluidas en el modelo.
RSE (Error estándar de los residuos): 0.7071. Este valor indica el tamaño promedio de los errores en las predicciones del modelo, siendo relativamente bajo, lo que sugiere que las predicciones son bastante precisas.
F-estadístico 32.97 con p-valor < 2.2e-16. Este resultado indica que el modelo es altamente significativo en su conjunto, lo que significa que al menos una de las variables predictoras tiene un impacto significativo en los niveles de PSA.
El modelo de regresión muestra que lcavol,
lweight, y svi son los factores más relevantes
para predecir los niveles de PSA, con relaciones significativas y
positivas. En cambio, lbph y lcp no tienen un
impacto significativo en los niveles de PSA. El \(R^{2}\) ajustado de 0.6248 indica que el
modelo explica el 62.48% de la variabilidad en PSA, y el RSE bajo
sugiere predicciones precisas. El modelo es útil, pero podría mejorarse
considerando otro conjunto de variables.
Para realizar las predicciones usando los modelos de
Ridge y Lasso utilizaremos la librería
glmnet. Primero preparamos nuestro conjunto de datos lo
dividimos entre el conjunto de training y el de
testeo. Para crear nuestro training set utilizaremos el
70% de las entradas de nuestro conjunto de datos.
library(glmnet)
# Preparamos un dataset sin las variables que no usamos
datos_usados <- subset(prostate, select = -c(age, pgg45, gleason))
# Dividimos el dataset entre predictores y el valor a predecir
X <- as.matrix(datos_usados[, -which(names(datos_usados) == "lpsa")])
y <- datos_usados$lpsa
# Creamos los conjuntos de train y test, especificando la seed para poder reproducir la ejecución
# con los mismos sets
set.seed(123)
train_idx <- sample(1:nrow(X), size = 0.7 * nrow(X))
X.train <- X[train_idx,]
y.train <- y[train_idx]
X.test <- X[-train_idx,]
y.test <- y[-train_idx]
Para utilizar el modelo de Ridge con glmnet debemos
poner el parámetro \(\alpha = 0\).
Utilizaremos un grid con 200 valores de \(\lambda\) diferentes, desde \(10^{4}\) hasta \(10^{-4}\). Estos valores se han determinado
después de realizar la gráfica varias veces para asegurarnos de que las
regiones que en ella se observan tienen una representación similar. A
partir del modelo, podemos obtener el conjunto de coeficientes y los
lambdas correspondientes a cada uno de los 200 modelos obtenidos.
lambdas = 10^seq(4,-4, length = 200)
ridge.model <- glmnet(X.train, y.train, alpha = 0, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
ridge.coefs <- as.matrix(ridge.model$beta)
ridge.lambda <- ridge.model$lambda
matplot(log(ridge.lambda), t(ridge.coefs), type = "l", lty = 1, col = rainbow(nrow(ridge.coefs)),
xlab = "log(λ)", ylab = "Coeficientes", main = "Ridge: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(ridge.coefs)), lty = 1, cex = 0.6)
Vemos claramente 3 regiones en la gráfica. Una primera región a la
izquierda, hasta \(\log\lambda \approx
-4\), donde los coeficientes son constantes. Una segunda región
intermedia, aproximadamente entre \(-4\) y \(4\) donde estos valores varían,
disminuyendo todos los coeficientes a excepción de lcp, el
cual comienza con un valor negativo y crece hasta tener uno positivo.
Finaliza con una tercera región, donde los coeficientes convergen a un
valor cercano a 0.
Obtenemos el valor óptimo estimado de \(\lambda\) usando validación cruzada, con la
función cv.glmnet con el parámetro \(\alpha=0\):
# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos
# especificar la seed para asegurarnos de que podemos reprocir la prueba
set.seed(123)
# Obtenemos el mejor valor de lambda usando cross-validation
cv.ridge.model <- cv.glmnet(X.train, y.train, alpha=0, lambda=lambdas)
lambda.min.ridge <- cv.ridge.model$lambda.min
El valor óptimo del \(\log\lambda_\text{Ridge}\) nos permite ubicar el valor del \(\lambda\) óptimo en la gráfica anterior, ubicándola en la región intermedia descrita con anterioridad.
A continuación presentamos la gráfica del \(MSE\) frente a \(\log\lambda\) correspondiente al modelo de Ridge.
plot(cv.ridge.model)
Vamos ahora a hacer la predicción de lpsa usando el
conjunto de testeo.
pred.ridge <- predict(cv.ridge.model, s = lambda.min.ridge, newx = X.test)
Para utilizar el modelo de Lasso con glmnet debemos
poner el parámetro \(\alpha = 1\).
Utilizaremos el mismo grid de usado para el modelo de
Ridge, pero con valores de \(\lambda\) entre \(10^1\) y \(10^{-4}\), valores que obtenemos de igual
manera tras reproducir la gráfica varias veces. Después procedemos de
forma similar.
lambdas = 10^seq(1,-4, length = 200)
lasso.model <- glmnet(X.train, y.train, alpha = 1, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
lasso.coefs <- as.matrix(lasso.model$beta)
lasso.lambda <- lasso.model$lambda
matplot(log(lasso.lambda), t(lasso.coefs), type = "l", lty = 1, col = rainbow(nrow(lasso.coefs)),
xlab = "log(λ)", ylab = "Coeficientes", main = "Lasso: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs)), lty = 1, cex = 0.6)
Podemos observar un gráfico similar al obtenido para Ridge,
con la diferencia de que Lasso sí elimina predictores, es decir
hace sus coeficientes 0. Contando esta diferencia, vemos también una
representación análoga de las 3 regiones que vimos en el gráfico del
modelo de Ridge, con una primera región donde el coeficiente de
lcp es negativo, hasta que se elimina, una región
intermedia donde se regulan el resto de parámetros, y una final donde se
eliminan todos los coeficientes (en Ridge también tendían a 0).
Es por esto que podemos suponer que, también en este caso, el valor
óptimo de \(\lambda\) se encontrará en
la región intermedia.
También podemos visualizar el plot por defecto del
modelo de Lasso, el cual nos muestra la evolución de los
coeficientes con el L1 Norm, que representa el error que se
intenta minimizar al utilizar este modelo:
\[ L1_\text{Norm} = RSS + \lambda \sum ^p_{j=1}|\beta_j| \]
plot(lasso.model, col = rainbow(nrow(lasso.coefs)))
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs),start=0), lty = 1, cex = 0.6)
Obtenemos el valor óptimo estimado de \(\lambda\) usando validación cruzada, con la
función cv.glmnet con el parámetro \(\alpha=1\):
# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos especificar la
# seed para asegurarnos de que podemos reprocir la prueba
set.seed(123)
# Obtenemos el mejor valor de lambda usando cross-validation
cv.lasso.model <- cv.glmnet(X.train, y.train, alpha=1, lambda=lambdas)
lambda.min.lasso <- cv.lasso.model$lambda.min
log(lambda.min.lasso)
## [1] -2.267873
A continuación presentamos la gráfica del \(MSE\) frente a \(\log\lambda\) correspondiente al modelo de Lasso.
plot(cv.lasso.model)
Vamos ahora a hacer la predicción de lpsa usando el
conjunto de testeo.
pred.lasso <- predict(cv.lasso.model, s = lambda.min.lasso, newx = X.test)
Creamos un data.frame que contenga los coeficientes para
cada modelo. Esta tabla nos sirve tanto para poder comparar
cuantitativamente los valores de cada coeficiente, como para preparar
los datos para presentarlos en un gráfico usando
ggplot.
coef.multi <- coef(modelo_lpsa)[-1]
coef.ridge <- coef(cv.ridge.model, s=lambda.min.ridge)[-1]
coef.lasso <- coef(cv.lasso.model, s=lambda.min.lasso)[-1]
coef.data <- data.frame(
variable = names(coef.multi),
multiple = as.numeric(coef.multi),
ridge = as.numeric(coef.ridge),
lasso = as.numeric(coef.lasso)
)
coef.data
## variable multiple ridge lasso
## 1 lcavol 0.54789698 0.33723614 0.3963463
## 2 lweight 0.53035882 0.55018591 0.4842218
## 3 lbph 0.07999311 0.04808135 0.0000000
## 4 lcp -0.03638017 0.57440441 0.5478204
## 5 svi 0.75974935 0.06876884 0.0000000
Cambiamos la tabla de datos a formato wide para poder
mapearlo correctamente usando ggplot
coef.data.long <- reshape2::melt(coef.data, id.vars = "variable",
variable.name = "Modelo", value.name = "Coeficiente")
ggplot(coef.data.long, aes(x = variable, y = Coeficiente, fill = Modelo)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Comparación de Coeficientes entre Modelos de Regresión",
x = "Predictor", y = "Peso del Coeficiente") +
theme_minimal()
lbph: Lo primero que salta a la
vista es que el modelo Lasso ha eliminado este
predictor a diferencia del resto, ya que el modelo
Lasso es el único que realiza selección de variables.
De forma similar, el resto de modelos le otorgan un peso muy
bajo.
lcavol: Los tres modelos coinciden
en otorgar un gran peso a este predictor, destacando el modelo lineal.
Lasso y Ridge otorgan un valor similar entre ellos, siendo el peso en el
modelo de Lasso mayor.
lcp: Este caso es especialmente
llamativo, y es que mientras los modelos Lasso y Ridge coinciden en
otorgarle un peso alto (es el coeficiente de mayor valor en ambos), en
el modelo de regresión múltiple este valor no es solo mucho menor, sino
que además es negativo.
lweight: En este predictor
encontramos el mayor consenso entre los tres modelos, otorgándole un
gran peso los tres modelos.
svi: Este es otro caso de gran
contraste entre los modelos, y es que mientras es el valor con mayor
peso en el modelo múltiple, con una gran distancia al segundo, el modelo
de Ridge le otorga un valor mínimo y el modelo de Lasso directamente lo
elimina.
El comportamiento tan diferente de los coeficientes lcp
y svi podría indicar la existencia de una correlación de
importancia notable entre dichas variables.
Vamos a calcular la predicción del modelo lineal (el resto ya han sido calculados), y a continuación calcularemos y compararemos el MSE (error cuadrático medio) que resulta de cada modelo.
# Funciones para calcular las métricas
mse <- function(y,y2) {
se = (y-y2)^2
mean(se)
}
rsq <- function(y,y2) {
rss = sum((y-y2)^2)
tss = sum((mean(y)-y)^2)
1 - rss / tss
}
rsq.adj <- function(y,y2,p) {
n = length(y)
1 - ((1-rsq(y,y2))*(n-1))/(n-p-1)
}
pred.multiple <- predict(modelo_lpsa, newdata = as.data.frame(X.test))
resultados <- data.frame(
Metric = c("MSE", "R2", "R2adj"),
Multiple = c(
mse(y.test, pred.multiple),
rsq(y.test, pred.multiple),
rsq.adj(y.test, pred.multiple, 5)
),
Ridge = c(
mse(y.test, pred.ridge),
rsq(y.test, pred.ridge),
rsq.adj(y.test, pred.ridge, 5)
),
Lasso = c(
mse(y.test, pred.lasso),
rsq(y.test, pred.lasso),
rsq.adj(y.test, pred.lasso, 4) # p = 4 aquí ya que, como recordamos, lasso elimina lbph
)
)
resultados
## Metric Multiple Ridge Lasso
## 1 MSE 0.4894444 0.6175785 0.6350265
## 2 R2 0.7134772 0.6384670 0.6282529
## 3 R2adj 0.6537850 0.5631476 0.5687733
Replicamos lo hecho para los coeficientes para graficar los resultados obtenidos:
resultados.long <- reshape2::melt(resultados, id.vars = "Metric", variable.name = "Modelo", value.name = "Valor")
ggplot(resultados.long, aes(Metric, Valor, fill = Modelo)) +
geom_bar(stat= "identity", position="dodge") +
labs(title = "Comparación de Métricas de Rendimiento entre Modelos",
x = "Métrica", y = "Valor de la Métrica") +
theme_minimal()
Como podemos observar, el modelo que mejor ha funcionado en nuestro caso es el de la regresión múltiple, teniendo un resultado mejor que el de los otros modelos tras medir su rendimiento usando MSE, \(R^2\) y \(R^2_\text{adj}\) en el conjunto de testeo.
Así, podemos concluir que los modelos que penalizan un número mayor de predictores, en este caso, funciona peor que el modelo de regresión múltiple que no cuenta con esa penalización. Además, tampoco representan una mejora en el rendimiento o en el espacio de almacenamiento, dado que contamos con un número muy limitado de predictores y de observaciones, por lo que no se recomendaría su uso en situaciones similares.
Para responder a si los modelos lineales podrían ser utilizados o no, revisamos los valores de las métricas obtenidas para el modelo de regresión múltiple:
Como vemos, nuestro modelo funciona relativamente bien con el conjunto de testeo, por lo que podemos concluir que este modelo podría llegar a ser utilizado para realizar predicciones.
suppressPackageStartupMessages({
library(MASS)
})
Ajustar el modelo LDA utilizando svi como variable de
respuesta y lcavol, lcp, lpsa
como variables predictoras.
prostate$svi <- as.factor(prostate$svi)
lda_model <- lda(svi ~ lcavol + lcp + lpsa, data = prostate)
El modelo muestra las probabilidades de cada grupo (0 y 1). Estas indican la proporción de datos que pertenece a cada clase.
Las medias de grupo presentan los valores medios
de las variables predictoras (lcavol, lcp,
lpsa) para cada clase del objetivo
svi.
Los coeficientes de los discriminantes lineales reflejan la importancia de cada variable en la construcción del discriminante lineal.
lda_model
## Call:
## lda(svi ~ lcavol + lcp + lpsa, data = prostate)
##
## Prior probabilities of groups:
## 1 2
## 0.7835052 0.2164948
##
## Group means:
## lcavol lcp lpsa
## 1 1.017892 -0.6715458 2.136592
## 2 2.551959 1.6018579 3.715360
##
## Coefficients of linear discriminants:
## LD1
## lcavol -0.08659598
## lcp 0.76640382
## lpsa 0.53153340
class: contiene las clases predichas por el
modelo.
posterior: contiene las probabilidades posteriores
de cada clase para cada observación.
x: valores de los discriminantes lineales para cada
observación.
lda_model.pred <- predict(lda_model)
names(lda_model.pred)
## [1] "class" "posterior" "x"
Las primeras tres observaciones fueron clasificadas en la clase 0, lo
que indica que el modelo ha predicho la ausencia de svi
para estos casos.
lda_model.pred$class[1:3]
## [1] 1 1 1
## Levels: 1 2
La matriz de confusión muestra que el modelo clasificó correctamente 69 observaciones del grupo 0 y 17 del grupo 1.
Hay 7 casos del grupo 1 que fueron clasificados incorrectamente como 0, y 4 casos del grupo 0 clasificados incorrectamente como 1.
Esto podría estar relacionado con un desbalance en las clases, ya que el grupo 0 tiene una representación mucho mayor (aproximadamente el 78%) que el grupo 1, lo que puede hacer que el modelo esté más inclinado a predecir el grupo mayoritario.
Este desbalance puede afectar la capacidad del modelo para detectar correctamente los casos del grupo minoritario, lo que es crucial en escenarios donde los casos de la clase menos frecuente son más importantes, como en la detección de enfermedades.
library(plotly)
conf_matrix <- table(lda_model.pred$class, svi)
conf_matrix_df <- as.data.frame(as.table(conf_matrix))
heatmap <- plot_ly(
data = conf_matrix_df,
x = ~Var1,
y = ~svi,
z = ~Freq,
type = "heatmap",
colors = c("white", "blue"),
showscale = FALSE
) %>%
layout(
title = "Predicciones de SVI",
xaxis = list(title = "Predicción"),
yaxis = list(title = "Clase Real")
)
heatmap %>%
add_annotations(
x = conf_matrix_df$Var1,
y = conf_matrix_df$svi,
text = conf_matrix_df$Freq,
showarrow = FALSE,
font = list(size = 12, color = "black")
)
Las probabilidades posteriores indican cuán probable es que cada observación pertenezca a cada una de las clases (0 o 1).
En este caso, para las primeras tres observaciones, la probabilidad de pertenecer a la clase 0 es extremadamente alta (cerca de 1).
lda_model.pred$posterior[1:3, ]
## 1 2
## 1 0.9998211 0.0001789316
## 2 0.9997230 0.0002769794
## 3 0.9997500 0.0002500003
Los valores de los discriminantes lineales representan las puntuaciones de cada observación en la función discriminante.
Estos valores son útiles para entender cómo el modelo discrimina entre los grupos basándose en las variables predictoras.
lda_model.pred$x[1:3]
## [1] -2.304200 -2.125721 -2.167584
El modelo LDA ha mostrado un rendimiento desigual debido a un posible desbalanceo de clases, donde el grupo 0 es significativamente más grande que el grupo 1. Esto ha llevado a que el modelo clasifique correctamente 69 observaciones del grupo 0, pero solo 17 del grupo 1. Aunque el número de falsos positivos es relativamente bajo (4), los falsos negativos (7) indican que el modelo no está detectando adecuadamente todas las observaciones del grupo 1. Este comportamiento es común en escenarios con clases desbalanceadas.
Aunque el modelo no es completamente ineficaz, su capacidad para identificar el grupo minoritario es limitada. Se podrían explorar técnicas como el ajuste el balanceo de clases (por ejemplo, sobremuestreo o submuestreo) para mejorar el rendimiento en la clasificación del grupo 1 y obtener una predicción más equilibrada.
svi en función de lcavol,
lcp y lpsaVeamos primero las distribuciones del svi frente a los
predictores propuestos:
boxplot1 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcavol,
type = "box", name="lcavol", width = 900) %>%
layout(xaxis = list(title = "svi"), yaxis = list(title = "lcavol"))
boxplot2 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcp,
type = "box", name="lcp", width = 900) %>%
layout(xaxis = list(title = "svi"), yaxis = list(title = "lcp"))
boxplot3 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lpsa,
type = "box", name="lpsa", width = 900) %>%
layout(xaxis = list(title = "svi"), yaxis = list(title = "lpsa"))
fig <- subplot(boxplot1, boxplot2, boxplot3, nrows = 1, titleX = TRUE, titleY = TRUE) %>%
layout(title="Distribución de los predictores según svi")
fig
Como podemos observar, la distribución de svi respecto
al resto de variables es muy diferenciable. A simple vista podemos ver
una clara relación entre valores negativos del svi con los
tres predictores propuestos.
Vamos ahora a calcular el modelo. Dado que la variable
svi es binaria (1 o 0), especificamos el parámetro
family=binomial:
# Nos aseguramos de que la variable svi esté en formato factor, por si se ha cambiado en algun apartado anterior
prostate$svi <- as.factor(prostate$svi)
modelo_logistico <- glm(svi ~ lcavol + lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico$coefficients, function(x) round(x,2))
El modelo resulta en la siguiente función:
\[ \text{svi} = -8.58 + 0.12 \cdot \text{lcavol} + 1.31 \cdot \text{lcp}+ 2.14 \cdot \text{lpsa} \]
Veamos el summary del modelo:
summary(modelo_logistico)
##
## Call:
## glm(formula = svi ~ lcavol + lcp + lpsa, family = binomial, data = prostate)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.5845 2.4942 -3.442 0.000578 ***
## lcavol 0.1209 0.7146 0.169 0.865645
## lcp 1.3113 0.4240 3.093 0.001984 **
## lpsa 2.1431 0.8037 2.667 0.007662 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 101.353 on 96 degrees of freedom
## Residual deviance: 39.624 on 93 degrees of freedom
## AIC: 47.624
##
## Number of Fisher Scoring iterations: 7
Lo primero que salta a la vista es que el error estimado para la
variable lcavol es unas 6 veces mayor que el coeficiente en
sí, lo cual nos indica que podría ser que la significancia de este
predictor sea baja. Además, recordemos del apartado 5
que lcavol guardaba una alta correlación con las otras
variables, lcp y lpsa. Esto se nos confirma al
observar el \(p-\)valor de esta
variable, de aproximadamente \(0.87\),
lo cual nos indica que este predictor no es estadísticamente
significativo.
Vamos a realizar una predicción de nuestro set original usando el modelo:
threshold = 0.5
log.probs = predict(modelo_logistico, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1
Veamos la tabla de valores originales:
table(log.preds)
## log.preds
## 0 1
## 79 18
Y a continuación la tabla de valores obtenidos:
table(prostate$svi)
##
## 1 2
## 76 21
Como se puede observar, el modelo solo falló en 3 ocasiones, dando un
\(0\) cuando debería ser \(1\), lo cual indica que el modelo hizo un
buen trabajo al predecir el conjunto de entrenamiento, con un 97% de
precisión. Es importante hacer notar que el dataset está cargado hacia
svi = 0, es decir, hay muchas más observaciones en las que
svi = 1, lo cual en general no es bueno para la generación
de modelos, y en particular puede ser que haya influenciado a esta
tendencia al 0 que muestra nuestro modelo.
De cualquier forma, en esta ocasión no hemos dividido el dataset en training y testing, por lo que estamos examinando el rendimiento de nuestro modelo en el mismo conjunto de datos con el que lo hemos entrenado. En este caso habría que tener en cuenta dos cosas:
Por estos motivos, se concluye que, a priori, no se puede decir que sea un buen modelo sin realizar primero pruebas de su rendimiento con un set de testeo, además de que mejoraría en gran medida si se entrenara con un set más equilibrado.
Vistos los resultados obtenidos con el modelo anterior, nos
preguntamos si podríamos mejorarlo más. Y es que, como explicamos, el
predictor lcavol tiene una baja significancia, esto junto
con la alta correlación con el resto de predictores del modelo nos lleva
a pensar que podemos mejorar el modelo al eliminar este predictor.
Vamos a proceder de forma similar al modelo original:
modelo_logistico.nolcavol <- glm(svi ~ lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico.nolcavol$coefficients, function(x) round(x,2))
summary(modelo_logistico.nolcavol)
##
## Call:
## glm(formula = svi ~ lcp + lpsa, family = binomial, data = prostate)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.5322 2.4718 -3.452 0.000557 ***
## lcp 1.3434 0.3830 3.507 0.000452 ***
## lpsa 2.1982 0.7416 2.964 0.003035 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 101.353 on 96 degrees of freedom
## Residual deviance: 39.653 on 94 degrees of freedom
## AIC: 45.653
##
## Number of Fisher Scoring iterations: 7
El modelo resulta en la siguiente función: \[ \text{svi} = -8.53 + 1.34 \cdot \text{lcp}+ 2.2 \cdot \text{lpsa} \]
Además, los \(p-\)values nos
muestran que todas las variables utilizadas, lcp y
lpsa en este caso, son significantes. Vamos a comparar los
resultados obtenidos
Vamos a realizar una predicción de nuestro set original usando el modelo:
threshold = 0.5
log.probs = predict(modelo_logistico.nolcavol, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1
Veamos la tabla de valores originales:
table(log.preds)
## log.preds
## 0 1
## 78 19
Y a continuación la tabla de valores obtenidos:
table(prostate$svi)
##
## 1 2
## 76 21
Observamos como hemos conseguido predecir un valor más utilizando
este modelo, que se desprende de la variable lcavol,
llegando así al 98% de precisión. Sin embargo, como también mencionamos
con el modelo anterior, debemos ser críticos con este modelo y prestar
especialmente atención al posible overfitting y al
desbalance del conjunto de datos.
Vamos a utilizar nuestro modelo para predecir un valor de
svi para valores concretos de los predictores
lcavol, lcp y lpsa. Utilizamos
los siguientes valores:
lcavol = 2.8269,lcp = 1.843,lpsa = 3.285Primero creamos el data.frame que contenga estos
datos:
pred_data = data.frame(
lcavol = 2.8269,
lcp = 1.843,
lpsa = 3.285
)
Realizamos la predicción usando nuestro modelo:
prob <- predict(modelo_logistico, newdata = pred_data, type = "response")
Obtenemos una probabilidad del 77.1% de que svi\(=1\) con estos valores de los
predictores.
Vamos a realizar la misma predicción, pero esta vez utilizando el
modelo propuesto en el subapartado anterior, que se desprendía del
predictor lcavol por su baja significancia estadística.
Dado que los datos ya están preparados, podemos calcular la probabilidad
con un solo comando (se ignora el valor de lcavol al
realizar la predicción):
prob.nolcavol <- predict(modelo_logistico.nolcavol, newdata = pred_data, type = "response")
Obtenemos una probabilidad del 76.2% de que svi\(=1\) con estos valores de los predictores.
Podemos observar que la probabilidad es prácticamente la misma que con
el modelo anterior.
summary del modeloPrimero, creamos nuestro modelo de componentes principales PCA:
# Creamos el subset para el PCA
pca.data <- prostate[,c("lcavol", "lweight", "lbph", "lcp")]
# Añadimos el svi a parte, ya que al ser categorico hay que reconvertirlo a numerico
pca.data["svi"] <- as.numeric(prostate$svi)
# Realizamos el PCA
pca.model <- prcomp(pca.data, scale. = TRUE, center = TRUE)
pca.var <- pca.model$sdev ^2
pca.pve <- pca.var / sum(pca.var)
summary del modelosummary(pca.model)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.5341 1.1862 0.7258 0.66771 0.51651
## Proportion of Variance 0.4707 0.2814 0.1054 0.08917 0.05336
## Cumulative Proportion 0.4707 0.7521 0.8575 0.94664 1.00000
Podemos observar en la primera fila del summary de
nuestro modelo, Standard deviation, la cantidad de variabilidad
que cada componente principal captura de los datos. Como se puede
observar, los valores disminuyen con cada componente, siendo la mayor
PC1 y la menor PC5.
La desviación estándar en un análisis de componentes principales representa la cantidad de variabilidad de los datos que es capturada por esta componente, lo cual es una medida clave en el PCA, ya que su objetivo es el de identificar las direcciones, o componentes principales, que mejor capturen la varianza del conjunto de datos.
El primer componente principal PC1 es, por definición,
aquel que caputra la mayor varianza posible, maximizando la varianza en
las observaciones proyectadas en esa dirección, con lo cual una gran
parte de las diferencias entre observaciones se ven reflejadas a lo
largo de esta dirección.
Por tanto, se puede decir que la mayor desviación estándar de esta primera componente es consecuencia del objetivo de un PCA.
Las siguientes filas del summary ahondan en esta
descripción, proporcionándonos valores porcentuales y acumulativos de la
descripción de la varianza explicada por cada componente. Representemos
estos valores gráficamente:
par(mfrow = c(1, 2))
plot(pca.pve, xlab = "Principal Component",
ylab = "Proportion of Variance Explained", ylim = c(0, 1),
type = "b")
plot(cumsum(pca.pve), xlab = "Principal Component",
ylab = "Cumulative Proportion of Variance Explained",
ylim = c(0, 1), type = "b")
Como podemos observar, un \(75%\) de la varianza de los datos originales se podrían explicar usando tan solo las dos primeras componentes principales, lo cual sugiere que estas dos componentes podrían ser suficientes para representar los datos, con una pequeña pérdida de información sobre la variación de los datos, llegando al \(85%\) en caso de utilizar tres componentes.
biplot con las dos primeras componentes
principalesRealicemos un biplot para ver cómo se expresarían nuestros datos en función de esas dos componentes.
biplot(pca.model, scale = 0)
En este biplot podemos observar los puntos que representan a los
pacientes en el espacio definido por las dos primeras componentes
principales (PC1 y PC2).
Podemos observar diferentes agrupaciones de puntos, lo que indicaría que estas observaciones tendrían valores similares también en las variables originales.
Además, gracias a las flechas que indican la composición de las
componentes principales según las variables originales, podemos observar
que, al haber un ángulo pequeño entre ellas, las variables
lcavol, lcp y svi están
fuertemente correlacionadas, como se corrobora en apartados anteriores.
Lo mismo ocurriría con lbph y lweight aunque
en menor medida, al haber un ángulo de separación mayor entre ellas.
Iniciamos con una carga de las librerias necesarias
library(caret)
Iniciamos echando una pequeña vista a la correlación existente entre los datos. Se ha eliminado la diagonal para mejorar la visualización.
df <- prostate
df$svi <- as.numeric(df$svi)
corrplot(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]),
type = "upper",
tl.cex = 0.7,
diag = FALSE,
addCoef.col = "black")
Si hubiéramos usado la función findCorrelation, no se
vería ninguna hasta modificar el valor por defecto de
cutoff = 0.95.
findCorrelation(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]),
names = TRUE,
cutoff = 0.44)
## [1] "lcavol" "lcp" "lweight"
suppressPackageStartupMessages({
library(pls)
})
Para abordar la colinealidad, se ajusta un modelo de PCR y uso validación cruzada para evaluar el error de predicción.
Se incluye scale = TRUE para normalizar los predictores
y validation = "CV" para usar validación cruzada en la
evaluación:
set.seed(42)
pcr.fit = pcr(lpsa ~ lcavol + lweight + lbph + lcp + svi,
data = df,
scale = TRUE,
validation = "CV")
Para decidir el número óptimo de componentes principales, se traza
las puntuaciones de la validación cruzada utilizando
validationplot(). Visualmente ya se ve que el menor error
de validación cruzada es con M = 5.
Además se ve que el aumento de componentes favorece en los resultados. Haciendo uso completo de los pedidos en el enunciado. Aunque de forma relativa el uso de 1 a 4 componentes podría ser equivalente si quisiéramos usar menos.
validationplot(pcr.fit, val.type = "MSEP")
Fijándonos en el summary() podemos ver mismamente que en
M = 5 se captura el 100% de la varianza en los
predictores.
summary(pcr.fit)
## Data: X dimension: 97 5
## Y dimension: 97 1
## Fit method: svdpc
## Number of components considered: 5
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps
## CV 1.16 0.7694 0.7558 0.7551 0.7567 0.7281
## adjCV 1.16 0.7679 0.7547 0.7541 0.7548 0.7258
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps
## X 47.07 75.21 85.75 94.66 100.00
## lpsa 56.48 58.88 59.44 60.56 64.43
O a través de la función which.min() vemos el error más
pequeño:
which.min( MSEP( pcr.fit )$val[1,1, ] ) - 1
## 5 comps
## 5
# Predicciones
pcr.pred_full <- predict(pcr.fit, df, ncomp = 5)
multi_pred_full <- predict(modelo_lpsa, newdata = df)
# SSE para PCR
sse_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)
# SSE para Regresión Lineal Múltiple
sse_multi_full <- sum((multi_pred_full - df$lpsa)^2)
n <- nrow(df)
p_pcr <- 5
p_multi <- length(coef(modelo_lpsa))
# Calcular el RSE para PCR
rse_pcr_full <- sqrt(sse_pcr_full / (n - p_pcr))
# Calcular el RSE para Regresión Lineal Múltiple
rse_multi_full <- sqrt(sse_multi_full / (n - p_multi))
# R2 pcr
ss_total_full <- sum((df$lpsa - mean(df$lpsa))^2)
ss_residual_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)
r_squared_pcr_full <- 1 - (ss_residual_pcr_full / ss_total_full)
adj_r_squared_pcr_full <- 1 - ((1 - r_squared_pcr_full) * (nrow(df) - 1) / (nrow(df) - 2 - 1))
# R2 multi
adj_r_squared_multi_full <- summary(modelo_lpsa)$adj.r.squared
cat("PCR: regresión de componentes principales\n",
"MULTI: Regresión Lineal Múltiple\n\n",
"############ RSE ############\n\n",
"RSE (PCR):", rse_pcr_full, "\n",
"RSE (MULTI):", rse_multi_full, "\n\n",
"############ R² ############\n\n",
"R² ajustado (PCR):", adj_r_squared_pcr_full, "\n",
"R² ajustado (MULTI):", adj_r_squared_multi_full, "\n")
## PCR: regresión de componentes principales
## MULTI: Regresión Lineal Múltiple
##
## ############ RSE ############
##
## RSE (PCR): 0.7032094
## RSE (MULTI): 0.7070626
##
## ############ R² ############
##
## R² ajustado (PCR): 0.6367798
## R² ajustado (MULTI): 0.6248055
RSE (Error Estándar de la Regresión):
Ambos modelos, PCR y Regresión Lineal Múltiple (MULTI), tienen un RSE muy similar:
Esto indica que ambos modelos tienen un error estándar de la regresión comparable, lo que sugiere que la dispersión de los errores de las predicciones es casi igual en ambos casos.
R² ajustado:
Conclusión:
Aunque ambos modelos presentan RSE similares, lo que refleja una dispersión de errores comparable, el modelo PCR tiene una ligera ventaja en cuanto a R² ajustado, indicando que PCR explica mejor la variabilidad de los datos. Este pequeño beneficio en el ajuste hace que PCR sea el modelo preferido en este caso, ya que demuestra un mayor rendimiento al ajustar los datos.
Sin embargo, la diferencia no es tan significativa, por lo que dependiendo del contexto y los requisitos del análisis, Regresión Lineal Múltiple (MULTI) también puede ser una opción válida si se busca simplicidad o interpretabilidad adicional.